%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, displaySciPy
Disclaimer: this course is adapted from the notebooks by
Introduction
SciPy is a scientific library that builds upon NumPy. Among others, SciPy deals with:
- Integration (scipy.integrate)
- Optimization (scipy.optimize)
- Interpolation (scipy.interpolate)
- Fourier Transform (scipy.fftpack)
- Signal Processing (scipy.signal)
- Linear Algebra (scipy.linalg)
- Sparse matrices (scipy.sparse)
- Statistics (scipy.stats)
- Image processing (scipy.ndimage)
- IO (input/output) (scipy.io)
Animation display with python
Animation with matplotlib
You can use FuncAnimation to animate a sequence of images:
%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")
def init():
ax.set_xlim(0, 2 * np.pi)
ax.set_ylim(-1, 1)
return (ln,)
def update(frame):
xdata.append(frame)
ydata.append(np.sin(frame))
ln.set_data(xdata, ydata)
return (ln,)
ani = FuncAnimation(
fig,
update,
interval=50,
frames=np.linspace(0, 2 * np.pi, 100),
init_func=init,
blit=True,
)display(HTML(ani.to_jshtml()))Another way of displaying video exists, using html5 video:
display(HTML(ani.to_html5_video()))References:
- Matplotlib Animations / JavaScript Widgets by Louis Tiao
Animation with plotly
matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:
import plotly.express as px
from plotly.offline import plot
df = px.data.gapminder()
fig = px.scatter(
df,
x="gdpPercap",
y="lifeExp",
animation_frame="year",
animation_group="country",
size="pop",
color="continent",
hover_name="country",
log_x=True,
size_max=55,
range_x=[100, 100000],
range_y=[25, 90],
)
fig.show("notebook")Linear algebra
scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.
References:
Solving linear systems:
Find x such that: A x = b for specified matrix A and vector b.
A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
b = np.array([[1, 2, 3]], dtype=np.float64).T
print(A, b)
x = linalg.solve(A, b)
print(x, x.shape, b.shape)[[ 1. 0. 3.]
[ 4. 5. 12.]
[ 7. 8. 9.]] [[1.]
[2.]
[3.]]
[[ 0.8 ]
[-0.4 ]
[ 0.06666667]] (3, 1) (3, 1)
Check the result at a given precision (different from ==)
np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)True
Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!
Eigenvalues/ Eigenvectors
A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:
A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)
np.allclose(A, evecs @ np.diag(evals) @ evecs.T)[-0.94760646+0.j 2.04440479+0.j 3.47941641+0.j]
------
[[ 0.50787323 -0.78900444 -0.34574091]
[-0.85203941 -0.51920701 -0.06673021]
[ 0.12686067 -0.32847537 0.93595422]]
True
Symmetric matrices
If A is symmetric you should use eigvalsh (H for Hermitian) instead: This is more robust and leverages the structures (you know they are real!)
Matrix operations
linalg.trace(A)# tracelinalg.det(A)# determinantlinalg.inv(A)# Inverse, consider NEVER using it though :)
Norms
print(linalg.norm(A, ord="fro")) # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))4.145345283920557
4.145345283920557
3.4794164120004525
[0.94760646 2.04440479 3.47941641]
4.144258214937556
A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))3.016294053009224
Random generation, distributions, etc.
References:
- Good practices with numpy random number generators by Albert Thomas
- Numpy documentation on RandomState
- Random Widgets, by Joseph Salmon: Visualization of various popular distributions.
seed = 12345
rng = np.random.default_rng(seed) # can be called without a seed
rng.random()0.22733602246716966
Optimization
Goal: find functions minima or maxima
References:
- Scipy Lectures on mathematical optimization.
from scipy import optimizeFinding (local!) minima
def f(x):
return 4 * x ** 3 + (x - 2) ** 2 + x ** 4
def mf(x):
return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)
xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()
Default solver for minimization/maximization: fmin_bfgs (see Wikipedia on BFGS)
x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)
plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()Optimization terminated successfully.
Current function value: -3.506641
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8
Optimization terminated successfully.
Current function value: -6.201654
Iterations: 5
Function evaluations: 12
Gradient evaluations: 6
Optimization terminated successfully.
Current function value: 2.804988
Iterations: 7
Function evaluations: 16
Gradient evaluations: 8

Find the zeros of a function
Find x such that f(x) = 0, with fsolve.
omega_c = 3.0
def f(omega):
return np.tan(2 * np.pi * omega) - omega_c / omega
x = np.linspace(1e-8, 3.2, 1000)
y = f(x)
# Remove vertical lines when the function flips signs
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)
optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))
my_zeros = (
np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation
from scipy.optimize import curve_fit
def f(x, a, b, c):
"""f(x) = a exp(-bx) + c."""
return a * np.exp(-b * x) + c
x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5) # true signal
yn = y + 0.2 * np.random.randn(len(x)) # noisy added
plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()
(a, b, c), _ = curve_fit(f, x, yn)
print(a,"\n", b,"\n", c)
2.6616069003229623
1.4624690685428108
0.5538807148955575
Displaying
plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()
For polynomial fitting, one can directly use numpy functionsL
x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:
Polyfit may be poorly conditioned

Interpolation
from scipy.interpolate import interp1d, CubicSpline
def f(x):
return np.sin(x)
n = np.arange(0, 10)
x = np.linspace(0, 9, 100)
y_meas = f(n) + 0.1 * np.random.randn(len(n)) # add noise
y_real = f(x)
linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)
cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)
plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()
Images
from scipy import ndimage, datasets
img = datasets.face()
type(img), img.dtype, img.ndim, img.shape
print(2 ** 8) # uint8-> code sur 256 niveau.
n_1, n_2, n_3 = img.shape
np.unique(img)
# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()256

RGB decomposition
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))
plt.tight_layout()
print(img.flags) # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags) # can edit now
img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show() C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : False
WRITEABLE : False
ALIGNED : True
WRITEBACKIFCOPY : False
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False

Convert a color image in grayscale
plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()
Blurring (🇫🇷: floutage)
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)
#| eval: false
%matplotlib widget
import ipywidgets as widgets
# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()
@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
"""Remove old lines from plots and plot new ones."""
# [l.remove() for l in ax.lines]
img_flou = ndimage.gaussian_filter(img, sigma)
ax.imshow(img_flou, cmap=plt.cm.gray)Changing colors in an image
import pooch
import os
url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)
img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)
fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape
ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)
ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)
plt.tight_layout()

References:
